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A numerical  procedure  has  been  developed  to  compute  the 
transonic  viscous -inviscid  interacting  flow  field  about  airfoils 
with  leading  edge  slats  or  trailing  edge  flaps.  The  inviscid 
theory  utilizes  conformal  mappings  in  a full  potential  flow 
formulation,  with  analytic  removal  of  singularities  and  a mixed 
flow  relaxation  procedure.  The  turbulent  boundary  layer  is 
computed  with  methods  based  on  a turbulent  kinetic  energy  form- 
ulation and  is  coupled  to  the  inviscid  flow  using  surface 
source  flow  boundary  conditions.  Semiempirical  methods  are 
used  to  compute  local  strong  interaction  regions  occurring 
in  the  vicinity  of  shock  waves,  trailing  edges,  and  to  account 
for  flow  separation.  Our  results  are  in  good  agreement  with 
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1.  INTRODUCTION 


For  an  aircraft  to  maneuver  effectively  in  the  transonic 
speed  range,  its  wing  must  generate  high  lift  coefficients 
without  incurring  excessive  drag  or  buffet.  The  recent 
development  of  supercritical  wings  can  enable  a designer  to 
meet  this  specification.  However  this  requirement  will  often 
degrade  the  aircraft’s  performance  at  cruise  with  larger  than 
optimal  drag  coefficients.  The  implementation  of  high-lift 
devices  at  transonic  speeds  offers  the  possibility  of  greatly 
enhancing  the  maneuvering  capabilities  of  modern  aircraft 
without  compromising  their  cruising  efficiency.  This  possibility 
has  been  proven  in  the  last  few  years  by  the  installation  of 
slats  on  the  F-4  and  the  positive  test  of  a slatted  wing  on  the 
F-14  aircraft.  The  performance  of  these  aircraft  in  managing 
climbs  and  turns  at  transonic  speeds  was  remarkably  improved 
by  the  presence  of  the  slats,  even  though  these  configurations 
have  not  been  shown  to  be  optimal  by  any  means. 

The  aerodynamic  designer  currently  lacks  an  analytical 
tool  to  design,  or  even  analyze,  transonic  airfoils  with 
high- lift  devices.  Furthermore  the  paucity  of  experimental 
data  currently  makes  it  difficult  to  determine  what  can  be 
achieved  with  these  maneuvering  devices.  Also,  accumulating 
experimental  data  on  such  configurations  would  be  extremely 
expensive  in  light  of  the  number  of  configurations  that  need 
be  tested  and  the  high  speeds  and  Reynolds  numbers  required  in 
a meaningful  wind-tunnel  test  program.  A theoretical  tool 
for  the  analysis  of  the  transonic  flow  over  two-element  airfoil 
systems  would  be  a valuable  first  step  in  aiding  the 
designer  by  reducing  the  number  of  configurations  that  need  be 
tested  and  by  providing  insight  into  the  flow  phenomena  that 
are  present  at  high  speeds. 

This  report  describes  the  development  of  a method  for 
niomerically  computing  the  viscous  transonic  flow  over  an 
airfoil  with  a leading-edge  slat  or  a trailing-edge  flap. 
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Practical  designs  will  primarily  require  slat  configurations, 
but  since  the  method  is  applicable  to  general  two-element 
airfoil  systems,  both  slats  and  flaps  are  considered.  In 
general,  even  the  inviscid  flow  fields  on  these  configurations 
are  difficult  to  obtain  analytically  because  of  the  complicated 
geometry  of  the  multiply-connected  domain.  Small-disturbance 
approximations  (such  as  that  used  for  this  problem  in  Ref.  1) 
do  not  appear  to  be  adequate  since  the  interaction  of  the  flow 
between  the  airfoils  will  provide  large  perturbations  to  the 
flow  field. 

In  recent  years  the  application  of  mixed-flow  relaxation 
techniques,  introduced  by  Murman  and  Cole  (Ref.  2)  has  made 
possible  the  numerical  computation  of  inviscid  transonic  flows 
over  a variety  of  geometrical  shapes  in  both  two  and  three 
dimensions.  These  methods  are  generally  based  on  the 
assiamption  of  irrotational  flow  and  solve  either  the  full 
potential  equation  or  an  appropriate  form  of  the  small-dis- 
turbance equations.  For  two-dimensional  flows  in  particular, 
accurate  and  efficient  solutions  to  the  full  potential 
equation  have  been  obtained  for  transonic  flows  over  airfoil 
sections  (Refs.  3 and  4)  over  axisymmetric  bodies  (Refs.  5 and  6) 
and  over  nacelles  (Refs.  7 and  8). 

More  recently,  as  reported  in  Refs.  9 and  10,  these 
relaxation  techniques  have  been  applied  to  compute  the  flow 
about  an  airfoil  with  a slat  or  a flap  at  transonic  speeds. 

The  approach  as  discussed  in  Ref.  10  and  summarized  in  Section 
2 is  to  solve  the  full  inviscid,  irrotational  flow  equations 
about  two-element  airfoil  systems.  The  methodology  consists 
of  the  1)  development  of  a suitable  computational  plane  and 
grid  system,  2)  evaluation  of  an  appropriate  set  of  governing 
inviscid  equations  and  boundary  conditions  in  terms  of 
smoothly- varying,  s ingle-valued  functions  in  the  computational 
domain,  and  3)  establishment  of  a stable  and  accurate  numerical 
procedure  for  the  solution  of  the  governing  equations.  An 
abbreviated  version  of  the  inviscid  analysis  appears  in 
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Ref.  9 with  further  details,  particularly  on  the  mapping 
methods  in  Ref.  10.  Arlinger  (Ref.  11)  has  recently 
independently  developed  a similar  inviscid  analysis  of  this 
problem. 

At  high  Reynolds  numbers,  solutions  of  the  inviscid 
flow  equations  provide  a reasonable  estimate  of  the  lift  on  an 
airfoil,  provided  the  angle  of  attack  is  below  that  for 
maximum  lift.  However,  inviscid  theory  provides  no  information 
on  skin  friction,  form  drag,  or  on  the  maximum  lift  of  an  airfoil. 
These  important  characteristics  are  completely  dominated  by 
boundary  layer  growth  on  airfoil  surfaces.  For  standard  airfoils 
at  low  speeds,  boundary  layer  effects  are  relatively  weak  at 
high  Reynolds  number  and  can  be  treated  as  a small  correction 
to  the  inviscid  solution.  At  transonic  speeds,  for  super- 
critical airfoils  and  multielement  airfoil  systems,  the 
situation  is  much  more  severe,  with  viscous  effects  playing 
a significant  role  in  reducing  the  lift  from  inviscid  values. 

Our  approach  to  computing  viscous  effects  on  two-element 
airfoil  systems  is  based  on  second-order  boundary  layer  theory 
with  a surface  source  formulation  of  the  viscous  matching 
conditions.  The  inviscid  and  viscous  flows  are  solved 
simultaneously  in  a self-consistent  fashion,  by  iteration.  The 
development  of  the  boundary  layer  over  the  surfaces  of  the 
airfoils  is  driven  by  the  inviscid  flow  with  the  equivalent 
mass  flow  boundary  condition  on  the  airfoil  surfaces.  The 
effect  of  the  boundary  layer  on  the  inviscid  flow,  and  the 
circulation  in  particular,  will  be  felt  through  the  renewed 
computation  with  equivalent  surface  sources.  In  our  analysis 
we  assume  that  the  airfoil  elements  are  sufficiently  far  apart 
so  that  the  boundary  layers  do  not  merge  in  the  slot  region. 

To  account  for  strong  interaction  regions  near  the  airfoil 
trailing  edges  and  in  the  vicinity  of  shock  waves,  semiempirical 
smoothing  procedures  are  used.  A rational  analytic  approach 
to  the  trailing  edge  interaction  has  recently  been  proposed  by 
Melnik,  Chow  and  Mead  (Ref.  12),  and  will  be  implemented  into 
our  approach  in  the  future.  Most  available  data  on  two-element 
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airfoil  systems  at  transonic  speeds  indicate  substantial 
regions  of  flow  separation.  In  the  absence  of  a definitive 
theoretical  method  for  treating  turbulent  separated  flows, 
we  resort  to  a semiempirical  procedure  to  model  this  phenomena. 
Details  of  our  viscous  flow  method  are  given  in  Section  3, 
with  a discussion  of  the  coupling  of  the  inviscid  and  boundary 
layer  flows  in  Section  4. 

We  have  applied  our  computation  to  several  typical 
slatted  and  flapped  airfoil  configurations  and  have  compared 
our  results  with  existing  wind  tunnel  test  data.  A discussion 
of  these  results  along  with  our  future  recommendations  of 
further  research  in  this  area  are  presented  in  Section  5. 


2.  INVISCID  SOLUTION 


MAPPINGS 

A crucial  step  in  the  development  of  a finite-difference 
method  to  compute  flows  over  complicated  geometries  is  to 
develop  a suitable  grid  system.  It  is  highly  desirable  to  have 
the  geometric  contours  aligned  with  a coordinate  line  in  order 
to  avoid  interpolations  and  extrapolations  in  applying  the 
surface  boundary  conditions.  It  is  also  convenient  for  external 
flow  problems  to  map  the  infinite  physical  domain  to  a finite 
computational  space  in  order  to  accurately  apply  the  far-field 
boundary  conditions.  Furthermore,  the  mappings  should 
concentrate  grid  lines  in  regions  of  steep  flow  gradients  such 
as  airfoil  leading  and  trailing  edges  and  in  the  slot  formed 
between  the  main  airfoil  and  the  slat  or  flap. 

In  our  approach,  we  use  analytic  and  numerical  conformal 
mapping  methods  to  transform  the  infinite  domain  external  to 
two-element  airfoil  system  to  the  annular  region  between  two 
concentric  circular  rings.  The  outer  ring  corresponds  to  the 
main  airfoil  surface  and  the  inner  ring  to  the  secondary  airfoil 
surface  (flap  or  slat).  Infinity  in  the  physical  plane  is 
mapped  to  a single  point  within  the  circular  annulus  in  the 
computational  domain.  The  mapping  method  follows  from  the  work 
of  Ives  (Ref.  13)  and  utilizes  a sequence  of  five  conformal 
transformations,  three  analytic  and  two  numerical.  The  mapping 
proceeds  as  follows;  First  the  main  airfoil  is  transformed  to 
a near  circle  by  a Von  Karman-Tref f tz  transformation.  This  is 
followed  by  a Theodorsen  transformation  utilizing  fast-Fourier 
transforms  to  map  the  main  airfoil  near  circle  to  an  exact 
circle.  The  third  mapping,  as  outlined  in  Ref.  13,  is  an 
analytic  transformation  of  the  secondary  airfoil  to  a near 
circle,  which  keeps  the  image  of  the  main  airfoil  a circle 
(but  of  different  radius).  (In  the  application  of  this  mapping, 
we  have  developed  an  approach  which  greatly  simplifies  some 
procedures  in  Ref.  13,  and  these  are  discussed  in  Ref.  10.) 
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Next,  the  near  circular  image  of  the  secondary  airfoil  is 
centered  at  the  interior  of  the  circle  corresponding  to  the  main 
airfoil  through  a bilinear  transformation.  And  finally,  the 
two  concentric  shapes  are  mapped  to  two  circular  rings  through 
a second  application  of  the  Theodorsen  transformation. 

An  orthogonal  grid  system  is  produced  by  taking  a polar- 
coordinate  system  (r,6)  emanating  from  the  center  of  the  circular 
annulus.  The  surfaces  of  the  two  airfoils  are  obtained  as  two 
constant  radius  lines,  r = 1 for  the  main  airfoil  and  r - r 

s 

for  the  secondary  airfoil.  The  point  corresponding  to  infinity 
in  the  physical  domain  is  located  at  r = r^  and  0=0.  A further 
analytic  coordinate  stretching  X = X (0)  is  used  in  the  circum- 
ferential direction  to  produce  a suitable  grid  spacing  in  the 
physical  plane  with  mesh  points  concentrated  near  leading 
and  trailing  edges  and  with  each  trailing  edge  coinciding 
with  a grid  point.  The  circumferential  stretching  discussed  in 
Ref,  10  required  the  program  user  to  select  three  parameters  to 
achieve  the  desired  point  concentration.  An  alternate  stretching 
which  requires  no  inputs  from  the  user  is  given  in  the  Appendix. 

A radial  stretching  Y = Y(r)  is  used  to  locate  the  point  of 
infinity,  r = r^,  midway  between  the  circular  airfoil  rings. 

The  final  computational  domain  is  sketched  in  Fig.  1.  In 
this  plane  a uniform  grid  produces  the  mesh  distribution  in  the 
annular  domain  shown  in  Fig.  2 and  in  the  physical  domain  shown 
in  Fig.  3 for  a typical  slat  configuration.  The  mapping  produces 
a grid  where  each  airfoil  surface  is  a coordinate  line,  the 
trailing  edges  occur  on  mesh  points  and  with  a high  density  of 
grid  lines  in  the  slot  region  and  near  all  stagnation  points. 
Although  the  mapping  procedure  is  quite  complicated,  our 
computer  program  to  calculate  the  coefficients  of  all  the  terms 
generally  requires  less  than  10  seconds  on  an  IBM  370/168. 
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MAIN  AIRFOIL 

Fig.  1 Computational  Domain 

FORMULATION 

The  governing  equations  for  the  inviscid,  irrotational 
compressible  flow  about  the  two-element  airfoil  system  are 
written  in  the  computational  domain  using  the  metric  H,  of  the 
above-mentioned  mappings.  A potential  function  is  introduced 
into  these  equations.  Singularities  are  seen  to  arise  for 
several  reasons.  Firstly,  the  metric  of  the  mapping  becomes 
unbounded  at  infinity  (r=  r^,  9 = 0).  A study  of  the  mapping 
function  indicates  that  near  infinity 

dt  k 

^ - — s-  as  Z ^ r rn 


i d 

where  c = x + iy  in  the  physical  plane  and  Z = r e in 

the  annular  domain,  and  k is  a known  complex  constant  taken  to 
ik 

be  k = kj^  e 2.  The  metric  may  then  be  regularized  by 


(2) 


where 

f -2r^r  cos0  + xj-  (3) 

and  H is  a smooth  bounded  function  which  goes  to  unity  at 


Next,  the  potential  function  itself  becomes  unbounded  and 
multivalued  near  the  point  of  infinity  in  the  computational 
domain.  One  contribution  to  the  singular  potential  comes  from 
the  behavior  of  the  uniform  flow  near  infinity,  which  can  be 
shown  to  be  of  the  form  for  Z -*•  r_ 

(t)j^(r,e)  = Real 


= ^r  cos  (e-hx-k2)  - r^  cos  (a-k2)J 


(4) 


where  a is  the  angle  of  attack  of  the  free-stream  velocity  vector. 


A second  contribution  to  the  singular  potential  comes 
from  the  multivalued  nature  of  the  circulatory  flow  near  infinity 
in  the  annular  domain.  In  taking  a closed  circuit  about  infinity 
Z = r^,  the  potential  must  jump  by  2t\  times  the  circulation 
about  each  airfoil.  The  solution  for  the  circulatory  flow 
potential  valid  near  infinity  is  found  as  a solution  to  the 
Prandtl-Glauert  equations  (see  Ref.  14)  and  is  transformed  to 
the  computational  domain  as 


>2(r,0)  = - (rj^+r2)  tan 


-1 


X - tan  ^ 


(5) 


where 
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-a  + tr-tan 


-1 


/ r sin  0 \ 

V cos  B-rJ 


(6) 
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and  where  is  the  free- stream  Mach  number  and  and  T2  are  the 
circulation  constants  about  the  main  and  secondary  airfoils, 
respectively. 

To  obtain  a single-valued  reduced  potential  another  term, 
((i^.  must  be  introduced  so  that  any  closed  contours  about  indivi- 
dual airfoils  will  produce  the  required  circulation  jump.  This 
is  obtained  through  a term 

<^3  = - ^2^ 

A reduced  potential  function  G(r,9)  may  now  be  defined 
which  remains  bounded  and  single-valued  throughout  the  entire 
annular  domain  as 

G(r,e)  = - $2  -'I>3  (8) 


The  governing  equations  then  become 


and 


- 1+  mi-  (u^+v2) 


where 


A-W^ 


1-M  sin  e 

oo 


-M  sin2f? 

00 


1-M  sin  e 
00 


and  f and  B are  defined  in  Eqs.  (3)  and  (6),  respectively.  Also, 
a is  the  speed  of  sound,  y the  ratio  of  specific  heats,  and  v 
and  u the  radial  and  circumferential  velocity  components, 
respectively,  given  by 


Hkj  I ^ 


E.r  sin  6 + 

1 oo 


''0 


- (Gg  - Tj)  - (r^  + rj) 


E,  (r„  cos  e-  r)  + 


“0 


where 


r^  cos  (9+a-k2)  - 2r^r  cos  (a-k2) 


^1  r 2 

^2^  = ^ r cos  (9+a- 
+ r^cos (9-a+k2) J 
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-"i  r 

“i  = r-  y 


r sin(6+a-k2)-2r^r  sin  (a-k2) 


-r 


sin  (0-a+k2) 


(15) 


The  boundary  condition  of  the  vanishing  of  the  normal 
velocity  on  the  surface  of  each  airfoil  becomes 


V = 0 on  r = r. 


(16) 


and 


V = 0 on  r = 1 


(17) 


which,  from  Eq. 
surfaces  r = r 


(12)  specifies  the  normal  derivative  on  the 


and  r = 1 . 


With  the  definition  of  the  reduced  potential,  Eq.  (8),  it 
can  be  shown  that  the  solution  at  Z = r may  be  specified  as 


G(r  ,0)  = 0 


(18) 


And,  finally,  the  Kutta  condition  requires  the  vanishing  of  the 
tangential  velocity  uH  at  each  trailing  edge,  which,  from  Eq.  (13) 
gives  two  linear  equations  in  the  two  unknowns  and  r2 . 

A good  set  of  initial  conditions  for  the  reduced  potential 
can  be  obtained  from  the  incompressible  solution  for  the  flow 
over  two  circles  developed  by  Legally  (Ref.  15). 

FINITE-DIFFERENCE  PROCEDURE 
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The  numerical  formulation  of  the  Neumann  boundary-value 
problem  described  above  for  mixed  subsonic  and  supersonic  flow 
follows  from  the  ideas  and  techniques  developed  for  the  single 
airfoil  problem.  The  reduced  potential  equation,  Eq.  (9)  is 
solved  by  a successive  column  relaxation  algorithm  utilizing 
type-dependent  differencing  originated  by  Murman  and  Cole 
(Ref.  2).  Our  method  stems  from  techniques  developed  by 
Jameson  (Ref.  16).  Since  our  mapping  produces  a grid  system 
that  does  not  remain  aligned  with  the  streamwise  direction, 
it  is  necessary  to  use  a coordinate  invariant  or  "rotated" 
difference  scheme,  developed  by  Jameson.  Furthermore,  it  is 
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necessary  to  develop  sweep  directions  that  are  less  than  90°  to 
the  streamline  direction.  A suitable  set  of  sweep  directions, 
illustrated  in  Figs.  1 and  2,  consists  of  first  dividing  the 
circular  annulus  into  two  sections  divided  by  the  ring  r = r^ 
or  Y = %.  Circumferential  lines  interior  to  r = r surround  the 
secondary  airfoil,  and  circumferential  lines  exterior  to  r = r^ 
surround  the  main  airfoil.  Sweep  directions  are  from  the 
stagnation  point  to  the  trailing  edge  of  each  airfoil  in  the 
annulus . 

At  the  end  of  each  sweep  new  values  of  the  circulation 
are  computed  by  setting  uH  = 0,  and  therefore  setting  the 
bracketed  term  on  the  right  hand  side  of  Eq.  (13)  to  zero,  at 
each  trailing  edge  and  solving  the  two  resulting  equations  for 
and  T 2-  In  some  cases,  where  there  are  large  circumferential 
gradients  near  one  of  the  trailing  edges,  such  as  a shock  or 
expansion  fan  near  the  trailing  edge  of  a slat,  stability  has 
been  enhanced  through  underrelaxation  of  the  determination  of 
the  circulation  constants. 


3.  BOUNDARY  LAYER  CALCULATION 
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The  development  of  the  boundary  layer  over  the  surfaces 
of  the  airfoils  is  assumed  to  be  driven  by  the  inviscid  flow 
with  equivalent  surface  sources.  At  transonic  speeds  the 
growth  of  the  boundary  layers  on  the  upper  and  lower  surfaces 
of  an  airfoil  is  highly  unsyrametrical . High  aft-loadings  cause 
a rapid  thickening  in  the  airfoil  top  and  a thinning  on  the  bottom 
as  the  trailing  edge  is  approached.  The  net  effects  is  to 
produce  a strong  uncambering  of  the  "equivalent"  airfoil  shape 
which  leads  to  a sharp  reduction  in  lift.  This  uncambering 
effect  could  also  be  looked  at  as  a strong  upwash  at  the  rear 
of  the  airfoil.  The  effect  of  the  boundary  layer  on  the  inviscid 
flow,  and  the  circulation  in  particular,  will  be  felt  through 
the  renewed  computation  of  the  equivalent  source  strength  on  the 
body.  It  is  assumed  that  gap  sizes  are  large  enough  so  that 
the  boundary  layers  of  neighboring  surfaces  do  not  merge.  In 
addition  the  effect  of  a finite  thickness  wake  passing  over  the 
downstream  element  is  ignored.  In  light  of  the  practical 
sizes  of  slats  these  assumptions  are  not  significant. 

The  growth  of  the  boundary  layer  in  its  initial  laminar 
stage  is  computed  using  an  integral  method  based  on  the 
approach  of  Thwaites  (Ref.  17).  The  particular  formulation 
employed  is  that  described  by  Rott  and  Crabtree  (Ref.  18) 
who  by  the  use  of  the  Illingworth-Stewartson  transformation 
showed  how  the  compressible  laminar  flow  on  a surface  is  reduced 
to  an  equivalent  incompressible  flow  that  can  be  computed  by 
Thwaites'  original  method.  At  transonic  speeds  the  laminar  run 
on  an  airfoil  surface  is  usually  quite  short  and  it  was  felt  that 
sufficient  accuracy  would  be  obtained  with  an  integral  method. 
Transition  is  still  an  imperfectly  understood  phenomenon  and 
difficult  to  predict.  Several  empirical  criteria,  such  as 
Crabtree's  and  Michel's  (Ref.  19)  are  available.  Alternatively, 
the  point  of  transition  can  be  specified  in  the  computational 
method.  Since  transition  is  most  often  fixed  in  wind  tunnel 
test  this  feature  is  extremely  useful.  In  addition,  should 
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separation  be  predicted  while  the  boundary  layer  is  still 
laminar,  transition  is  again  assumed  to  occur,  and  the 
computation  is  continued.  Transition  is  assumed  to  occur 
instantaneously,  and  from  this  point  on  the  turbulent  boundary 
layer  solution  is  obtained  with  Bradshaw's  finite  difference 
method  (Ref.  20) . Starting  conditions  for  the  turbulent 
calculation  are  obtained  by  requiring  continuity  of  the  mass 
and  momentum  fluxes  within  the  boundary  layer  during  transiton. 

In  the  method  of  Bradshaw  three  equations,  the  mean  motion 
equations  for  continuity  and  momentum  and  an  empirical  equation 
for  the  shear  stress  obtained  from  the  exact  turbulent  energy 
equation,  are  integrated  numerically  along  the  surface.  Since 
the  three  equations  are  of  an  hyperbolic  nature  the  integration 
is  performed  by  marching  along  the  surface.  This  method  has  been 
shown  to  be  very  accurate  for  a wide  variety  of  flows.  Al- 
ternatively, the  turbulent  boundary  layer  growth  can  be  computed 
by  Green's  integral  lag-entrainment  method  (Ref.  21).  Green's 
method  solves  at  each  station  along  the  surface  a system  of  three 
equations,  the  momentum  integral  equation,  an  entrainment 
equation  and  an  equation  for  the  streamwise  rate  of  change  of 
the  entrainment  coefficient.  The  last  of  these  equations  was 
developed  from  Bradshaw' s empirical  equation  for  the  shear 
stress.  The  two  methods  have  thus  the  same  physical  foundations, 
and  the  results  of  the  two  methods  agree  very  well  and  produce 
essentially  the  same  results  in  the  program.  Both  methods  are 
included  in  tie  program  now  because  of  the  advaptages  each  might 
have  in  future  developments  of  the  viscous  analysis  method.  Green's 
method  is  capable  of  continuing  the  calculation  beyond  the 
trailing  edge  to  determine  the  thickness  of  the  wake.  At  a 
later  time  it  is  planned  to  examine  the  effects  of  a finite 
thickness  wake  and  of  the  merging  of  the  boundary  layers  and/or 
wakes  from  the  two  airfoil  elements  on  the  results  of  the  program 
especially  for  configurations  where  the  elements  are  closely 
spaced.  In  cases  where  the  interaction  of  the  merging  shear 
layers  is  strong  integral  methods  may  become  inaccurate.  Thus, 
even  though  at  present  Green's  method  is  sufficient,  the  possible 
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use  of  Bradshaw's  finite  difference  procedure  in  the  computation 
of  merging  flows  justifies  its  presence  in  the  program. 

The  parameters  of  interest  which  the  boundary  layer 
computation  yields  are  the  displacement  thickness  and  the  skin 
friction.  The  latter  is  integrated  to  give  the  skin  friction 
drag  on  the  airfoil  conf iguraton.  Following  the  approaches  taken 
in  the  theoretical  methods  of  analysis  for  single-airfoil  tran- 
sonic flows  and  multi-element  incompressible  flows,  the  major 
effect  of  the  boundary  layer  on  the  inviscid  flow  is  through  weak 
displacement  effects.  Within  these  theories  correction  to  the 
inviscid  solution  can  be  obtained  by  determining  the  inviscid 
flow  over  an  equivalent  body  obtained  either  by  adding  the  dis- 
placement thickness  to  the  airfoil,  or  alternatively,  as  in  the 
present  approach,  by  allowing  for  an  appropriate  mass  flow  at 
the  airfoil  surface  (Ref.  22). 

These  procedures  are  uniformly  valid  in  regions  where 
the  surface  geometry  is  smooth  and  the  inviscid  surface  pressures 


are  regular.  Unfortunately  in  regions  of  strong  interactions, 
such  as  trailing  edges  and  shocks,  ordinary  boundary  layer  theory 
breaks  down.  In  the  past  empirical  corrections  have  been  used 
quite  successfully  to  compute  the  displacement  effects  in  strong 
interaction  regions.  Since  the  purpose  of  constructing  an 
equivalent  body  is  to  have  a surface  around  which  the  flow  can 
be  considered  inviscid,  a way  to  account  for  strong  interaction 
regions,  short  of  actually  building  accurate  theoretical  models 
and  solving  them,  is  to  model  the  inviscid  streamline  in  these 
regions.  This  is  the  approach  that  is  followed  here. 

The  shock  wave-boundary  layer  interaction  has  somewhat  of 
a local  nature.  The  boundary  layer  spreads  the  pressure  gradient 
of  the  shock  in  the  free  stream  over  several  boundary  layer 
thicknesses  on  the  surface  (Ref.  23).  Often  a bubble  of  separated 
flow  will  occur  under  the  foot  of  the  shock.  In  such  a case, 
the  flow  around  the  airfoil  is  only  affected  locally.  Only  if 
the  shock  is  strong  enough  to  separate  the  boundary  layer  will 
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the  flow  at  the  trailing  edge  and,  hence,  the  circulation  be 
greatly  affected.  The  approach,  thus,  is  to  smear  the  pressure 
rise  through  the  shock  over  a short  distance  and  continue  the 
boundary  layer  computation  still  accepting  the  calculated  dis- 
placement thickness  as  an  adequate  representation  of  the  inviscid 
streamline.  Often  the  inviscid  calculation  itself  will  smooth 
the  shock  sufficiently  for  the  boundary  layer  computation  to 
proceed  without  separating.  A steep  pressure  rise  is  spread 
over  a short  distance  will  produce  a "bump"  in  the  displacement 
thickness.  As  the  rise  is  spread  over  a larger  distance  this 
"bump"  will  decrease  and  enough  smoothing  will  eliminate  it. 
Downstream  of  the  shock  wave  the  boundary  layer  characteristics 
are  remarkably  independent  of  the  degree  of  smoothing.  The 
downstream  boundary  layer  growth  is  practically  independent 
of  the  shock  modeling;  and  it  is  the  rate  of  growth  rather  than 
the  magnitude  of  the  displacement  thickness  as  the  trailing 


edge  is  approached  that  determines  the  decrease  in  circulation 
due  to  the  boundary  layer.  The  displacement  surface  obtained 
in  this  manner  near  the  shock  approximates  very  closely  the 
"bump"  model  proposed  by  Yoshihara  and  Murman  (private  commun- 
ication) without  a tendency  to  introduce  instabilities  in  the 
solution . 


Near  a trailing  edge  second-order  botindary  layer  theory 
again  becomes  singular.  The  displacement  thickness  slope  and 
hence  equivalent  mass  flow  grows  without  bound  as  the  trailing 
edge  is  approached.  Boundary  layer  calculation  programs  such 
as  Bradshaw's  or  Green's  reflect  this  feature.  There  is 
extensive  evidence  (Ref.  12)  to  indicate  that  the  trailing 
edge  singularity  is  eliminated  when  the  boundary  layer  and 
inviscid  flow  are  solved  simultaneously,  in  an  iterative  fashion. 
However  due  to  possible  numerical  errors  in  the  trailing  edge 
region,  the  displacement  thickness  is  smoothed  on  a scale  of  a 
few  boundary  layer  thicknesses  at  the  trailing  edge. 


The  approach  followed  in  single-element  airfoil  programs 
such  as  those  described  in  Refs.  24-26  has  been  to  extrapolate 
6*  from  a short  distance  upstream,  thus  providing  a smooth 
streamline  passing  over  the  trailing  edge.  This  method  gives  a 
good  representation  of  the  equivalent  body  locally.  As  such, 
this  is  the  approach  currently  followed  in  the  present  method. 
Recent  developments  will  make  available  an  approach  that  does 
not  rely  on  this  kind  of  empiricism.  The  work  of  Melnik,  Chow 
and  Mead  (Ref.  12)  provides  a closed  form  representation  for  the 
local  trailing  edge  solution,  and  this  would  replace  the 
extrapolation  procedure  currently  relied  on. 
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4.  COUPLING  OF  INVISCID  FLOW  AND  BOUNDARY  LAYER 


The  viscous  flow  over  the  two-element  system  is  computed 
in  the  form  of  an  inviscid  flow  over  an  equivalent  set  of 
surfaces.  These  are  the  streamlines  closest  to  the  airfoil 
elements  where  viscosity  can  be  neglected.  The  equivalent  body 
can  be  obtained  by  numerically  changing  the  ordinates  of  the 
airfoil  system,  which  would  necessitate  a remapping  of  the 
configuration  in  order  to  compute  a new  inviscid  flow.  An 
equivalent  and  simpler  procedure  which  avoids  the  need  to  remap 
is  to  change  the  boundary  condition  at  the  airfoil  surface  from 
one  of  zero  flow  through  the  surface  to  one  of  nonzero  normal 
velocity  with  a value  wich  makes  the  equivalent  body  a stream- 
line. In  this  formulation  a source  distribution  computed  from 
the  displacement  thickness  is  placed  at  the  surface  of  the 
airfoils  (i.e.  seeLighthill  (Ref.  22)).  This  source  distribution 
is  obtained  by  enforcing  the  condition  that  the  displacement 
surface  be  a streamline.  In  the  notation  of  Fig.  4 this 
means  that  the  streamfunction,  ij;  , defining  the  equivalent 
body  is  equal  to  zero.  The  nearby  airfoil  surface  (no  longer 
a streamline),  here  denoted  by  ipQ,  can  be  developed  from 
by  a one  term  Taylor  series  expansion;  thus 

'h  = ’^'e  • ^e  If  ^ ty  (19) 

Since  y = 6*  and  -1^  = +p.U  , it  follows  that 
e ay  e e’ 

'<'0  " ■ (20) 


Differentiating  with  respect  to  x yields 
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Now,  — — = -pV_,  where  V the  normal  velocity  component  at  the 

oX  0 ^ ^ 

airfoil  surface.  Thus, 

''s  = ? Is:  <Pe"e«*> 

The  new  boundary  conditions  to  be  satisfied  are  then 
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because  in  the  computational  plane  a positive  source  strength 
is  in  the  direction  of  negative  computational  velocities.  In 
terms  of  the  reduced  potential  these  boundary  conditions  become 
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The  iteration  process  begins  with  the  computation  of  the 
inviscid  flow  over  the  airfoil  systems.  The  first  calculation 
of  the  displacement  thickness  is  done  using  this  inviscid  flow. 
Successive  computations  are  made  using  the  inviscid  flow  over 
the  equivalent  bodies,  these  being  obtained  using  nonzero  Nevimann 
boundary  conditions  at  the  body  surfaces.  Computations  are 
done  iteratively  until  compatibility  between  inviscid  flow  and 
displacement  thickness  is  obtained,  usually  three  or  four 
iterations . 

Occasionally  in  the  first  computation- of  the  inviscid  flow, 
velocities  larger  than  the  limiting  velocity  are  generated  near  the 
leading  edge  of  the  slat.  To  give  the  iteration  cycle  a chance 
to  make  an  initial  correction  to  the  circulation  and  thus  cut 
down  on  these  expansions,  a lower  limit  on  the  speed  of  sound  is 
set.  This  allows  the  relaxation  to  proceed  with  finite  veloc- 
ities everywhere.  Also,  after  the  source  distribution  is  added 
on  the  surface  some  "wiggles"  may  develop  in  the  potential  flow 
solution  near  the  trailing  edge;  these  are  caused  by  the  large 
uncambering  of  the  airfoil  in  the  region  and  its  large  associated 
source  distribution  which  suddenly  ends  at  the  trailing  edge. 

Such  irregularities  in  the  pressure  distribution  would  cause 
similar  "wiggles"  in  the  subsequent  boundary  layer  computation 
and  one  could  not  expect  the  iteration  process  to  converge  under 
these  conditions.  Therefore  the  pressures  near  a trailing  edge 
are  extrapolated  from  a short  distance  upstream. 

Separation  is  unfortunately  a common  phenomenon  with  high- 
lift  configurations,  and  even  more  so  at  transonic  speeds. 

Presently  no  attempt  is  made  to  model  regions  of  separation  that 
do  not  envelope  the  trailing  edge;  such  separation  regions 
usually  occur  on  the  concave  lower  surface  of  a slat.  If  the  flow 
should  separate  before  the  trailing  edge  is  reached  and  the 
separation  region  passes  over  it,  an  attempt  to  model  the 
streamline  passing  over  the  separation  region  is  made  since  the 
direction  of  the  flow  near  the  trailing  edge  is  critical  to  the 
determination  of  the  circulation.  An  empirical  formula  relating 
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the  height  between  the  trailing  edge  and  the  "inviscid"  stream- 
line to  the  free  stream  conditions  on  local  geometry  is 
utilized  and  a linearly  growing  source  distribution  that  will 
generate  this  height  is  constructed. 
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5.  RESULTS  AND  DISCUSSION 
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The  method  has  been  applied  to  a variety  of  two-element 
airfoil  configurations  and  a few  typical  results  are  presented 
here.  In  order  to  evaluate  the  performance  of  the  method 
independent  of  the  particular  strong  interaction  models  and 
possible  wind  tunnel  blockage  corrections,  an  essentially 
incompressible  case  is  considered  first.  Figure  5 shows  the 
computed  and  measured  (Ref.  27)  surface  pressure  distributions 
on  a NACA  23012  airfoil  with  a 2H  flap.  The  plot  depicts  two 
inviscid  calculations,  the  viscous  interaction  calculation 
and  the  wind  tunnel  data.  The  computed  boundary  layer  growth 
on  the  main  airfoil  surface  is  small.  On  the  flap  upper  surface, 
a separated  flow  region  occurs  which  has  a large  effect  on  the 
lift.  The  agreement  with  experimental  data  is  excellent  except 
in  the  vicinity  of  the  leading  edge  of  the  flap.  The  discrepancies 
in  this  region  are  due  to  slight  differences  between  the  goemetry 
of  the  configuration  tested  and  that  modeled  by  the  computation. 

The  wind  tunnel  model  has  a long  lip  extending  from  the  trailing 
edge  of  the  main  airfoil.  This  protruberance,  whose  length  is 
about  57o  of  the  chord,  was  used  to  seal  off  the  slot  when  the 
flap  was  retracted  and  reached  over  the  leading  edge  of  the 
flap,  when  it  was  extended.  The  conformal  mapping  method  used  in 
the  computation  cannot  handle  this  geometric  complexity.  However, 
the  modeled  geometry,  as  shown  in  Fig.  5 seemed  to  produce  quite 
acceptable  results  over  most  of  the  configuration. 

Little  transonic  data  on  airfoils  with  leading  edge  slats 
is  available  and,  like  the  previous  case,  these  configurations 
have  regions  of  separated  flow.  In  Fig.  6,  the  computed  pressure 
distribution  and  the  experimental  data  (Ref.  28)  for  a NACA 
64A010  airfoil  with  a slat  at  = 0.7,  a = 6°  and  Re  = 7.8 
million  are  compared.  Lower  surface  separation  on  the  slat  alters 
drastically  the  flow  through  the  slot  and  again,  substantial 
discrepancies  occur  near  the  leading  edge  of  the  downstream 
element  (for  this  case  the  main  airfoil).  The  method  predicted 
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separation  near  the  lower  corner  of  the  slat  (which  has  been 
rounded  slightly)  but  no  attempt  was  made  to  model  the  massive 
separation  region  on  the  concave  surface  of  the  slat.  Separation 
was  also  predicted  on  the  upper  surface  of  the  slat,  near  the 
leading  edge.  On  the  main  airfoil,  away  from  the  leading  edge 
(slot)  region  satisfactory  agreement  with  the  data  is  obtained. 
Computed  constant  Mach  number  lines  are  shown  in  Fig.  7,  indi- 
cating the  imbedded  regions  of  supersonic  flow. 
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Fig.  7 Computed  Mach  Number  Contours  — NACA  64A010  Airfoil  with  18A 
Slat  - Moo  = 0-7.  a « 6® 

The  agreement  is  improved  when  the  slat  is  more  slender 
as  in  the  configuration  shown  in  Fig.  8.  The  arrangement  was 
obtained  from  a basic  NACA  64A406  airfoil  (Ref.  1).  The  flow 
separation  on  the  slat  does  not  affect  the  flow  coming  out  of 
the  slot  as  much  as  it  did  in  the  previous  case.  The  pressure 
distribution  on  the  main  airfoil  is  predicted  quite  well 
including  the  multiple  peaks  near  the  leading  edge.  Figure  9 
shows  the  Mach  number  contours  for  this  case.  It  is  interesting 
to  note  the  pocket  of  supersonic  flow  existing  in  the  slot. 

The  exit  of  the  slot  is  essentially  sonic  with  the  flow  quickly 
reaccelerating  to  supersonic  velocities  behind  it.  This  pattern 
is  reflected  in  the  multiple  peaks  in  both  the  computed  and 
experimental  pressures. 

Figures  10  and  11  give  results  for  two  different  positions 
of  the  slat.  Both  are  high  angle  of  attack  cases  and  stall 
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was  being  approached  in  both.  But  it  should  be  noted  that  at 
these  higher  incidences  the  amount  of  separation  on  the  lower 
surface  of  the  slat  is  substantially  reduced  and,  as  a result, 
the  agreement  between  data  and  computation  is  considerably 
improved  in  this  region. 


All  the  calculations  were  done  at  the  Mach  nvimber  and 
angle  of  attack  quoted  for  the  data.  Thus  the  possible 
influence  of  wind  tunnel  blockage  and  flow  angularities  was 
not  taken  into  accomt.  Not  unexpectedly,  the  comparisons  with 
data  have  been  marred  by  the  presence  of  substantial  regions 
of  separation.  But  the  reason  for  undertaking  this  study  was 
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APPENDIX 


CIRCUMFERENTIAL  STRETCHING 

An  alternate  circumferential  stretching  has  been  developed 
recently.  In  Ref.  10  it  was  reported  that  a suitable  grid  could 
be  generated  with  the  stretching 

e = C]^  sin  (J)  + C2  sin  2(|)  + c^  (A-1) 

where  <J)  is  the  azimuthal  location  of  a point  in  the  unstretched 
plane  and  9 is  the  stretched  coordinate,  Cj^,  C2,  and  c^  are 
three  constants  to  be  chosen  by  the  user  for  each  configuration. 

The  proper  concentration  of  grid  points  near  the  leading  and 
trailing  edges  can  be  obtained  also  by  having  the  azimuthal 
location  of  the  same  determined  by  the  intersection  of  the  ring  r = 1 
with  radial  lines  emanating  from  the  point  of  infinity  at  equal 
increments  as  shown  in  Fig.  A-1.  In  the  notation  of  Fig.  A-1, 
these  points  can  be  located  by 

cos  0 = r^  sin^t  + X + r^  sin^<|)  cos  ((>  (A-2) 

This  stretching  does  not  require  any  inputs  from  the  user. 

In  the  generation  of  the  grid  ({>  would  be  chosen  at  equal  increments. 
This  stretching  is  a natural  development  of  the  mapping  which 
transformed  infinity  to  a single  off-center  point.  For  the 
case  where  r and  r go  to  zero  9 and  d*  become  identical.  If  the 
point  of  infinity  were  to  be  moved  away  from  the  center  of  the 
outer  circle  the  same  mapping  would  pull  points  on  the  unit  circle 
around  the  circumference.  The  stretching  given  by  Eq.  (A-2) 
tends  to  reflect  this  phenomenon.  An  additional  stretching  which 
locates  the  point  of  infinity  and  each  trailing  edge  exactly 
on  grid  lines  is  still  used  and  has  been  described  in  Ref.  10. 
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